{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "private_outputs": true,
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/mikexcohen/Statistics_book/blob/main/stats_ch06_transformations.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Modern statistics: Intuition, Math, Python, R\n",
    "## Mike X Cohen (sincxpress.com)\n",
    "#### https://www.amazon.com/dp/B0CQRGWGLY\n",
    "#### Code for chapter 6 (transformations)\n",
    "\n",
    "---\n",
    "\n",
    "# About this code file:\n",
    "\n",
    "### This notebook will reproduce most of the figures in this chapter (some figures were made in Inkscape), and illustrate the statistical concepts explained in the text. The point of providing the code is not just for you to recreate the figures, but for you to modify, adapt, explore, and experiment with the code.\n",
    "\n",
    "### Solutions to all exercises are at the bottom of the notebook.\n",
    "\n",
    "#### This code was written in google-colab. The notebook may require some modifications if you use a different IDE."
   ],
   "metadata": {
    "id": "mzfXc9E3Xq7k"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dFT1qwVEyxTW"
   },
   "outputs": [],
   "source": [
    "# import libraries and define global settings\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# define global figure properties used for publication\n",
    "import matplotlib_inline.backend_inline\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats('svg') # display figures in vector format\n",
    "plt.rcParams.update({'font.size':14,             # font size\n",
    "                     'savefig.dpi':300,          # output resolution\n",
    "                     'axes.titlelocation':'left',# title location\n",
    "                     'axes.spines.right':False,  # remove axis bounding box\n",
    "                     'axes.spines.top':False,    # remove axis bounding box\n",
    "                     })"
   ]
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "IpajvgPyTI22"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Z-score"
   ],
   "metadata": {
    "id": "u5X-5wVwTI0P"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "X = [ 1,4,-5,2,1 ]\n",
    "\n",
    "mu = np.mean(X)\n",
    "sigma = np.std(X,ddof=1)\n",
    "\n",
    "X_z = (X-mu) / sigma\n",
    "print(X_z)"
   ],
   "metadata": {
    "id": "1EOjqcR0TIxp"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "print(f'       Original | z-transformed')\n",
    "print(f'Mean:      {np.mean(X):.2f} | {np.mean(X_z):.2f}')\n",
    "print(f'stdev:     {np.std(X,ddof=1):.2f} | {np.std(X_z,ddof=1):.2f}')"
   ],
   "metadata": {
    "id": "Rmpi9qVDTIvC"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "qWql-E8EdCKD"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.1: Example distributions of height and weight"
   ],
   "metadata": {
    "id": "PyHezTLsdCFW"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "### Note about this code: This code actually provides my solution to Exercise 6.\n",
    "# So if you want the full experience of the exercise, please don't inspect this code too carefully!"
   ],
   "metadata": {
    "id": "8hmBPJHqsFM4"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# fake heights and weights, in units of cm and kg\n",
    "N = 3425\n",
    "height = np.arctanh(np.random.uniform(-.9,.8,size=N)) * 20 + 160 + np.random.randn(N)*3\n",
    "weight = np.arctanh(np.random.uniform(-.3,.99,size=N)) * 10 +  70 + np.random.randn(N)*3\n",
    "\n",
    "# our imaginary individual\n",
    "ind_height = 177\n",
    "ind_weight = 70\n",
    "\n",
    "\n",
    "\n",
    "# z-score the distributions\n",
    "height_z = (height-np.mean(height)) / np.std(height,ddof=1)\n",
    "weight_z = (weight-np.mean(weight)) / np.std(weight,ddof=1)\n",
    "\n",
    "# z-score the individual\n",
    "ind_height_z = (ind_height-np.mean(height)) / np.std(height,ddof=1)\n",
    "ind_weight_z = (ind_weight-np.mean(weight)) / np.std(weight,ddof=1)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# figure layout\n",
    "_,axs = plt.subplots(2,3,figsize=(12,6))\n",
    "\n",
    "# plot the values\n",
    "axs[0,0].axvline(ind_height,color='k',linestyle='--',linewidth=2)\n",
    "axs[0,0].set(xlabel='Height (cm)',yticks=[],xlim=[ind_height-50,ind_height+50])\n",
    "axs[0,0].set_title(r'$\\bf{A}_1$)  Raw height')\n",
    "\n",
    "axs[1,0].axvline(ind_weight,color='k',linestyle='--',linewidth=2)\n",
    "axs[1,0].set(xlabel='Weight (kg)',yticks=[],xlim=[ind_weight-20,ind_weight+20])\n",
    "axs[1,0].set_title(r'$\\bf{A}_2$)  Raw weights')\n",
    "\n",
    "\n",
    "\n",
    "# plot the raw distributions with the individual\n",
    "axs[0,1].hist(height,bins='fd',color=(.5,.5,.5),edgecolor=(.8,.8,.8))\n",
    "axs[0,1].axvline(ind_height,color='k',linestyle='--',linewidth=2)\n",
    "axs[0,1].set(xlabel='Height (cm)',ylabel='Count',xlim=[ind_height-50,ind_height+50])\n",
    "axs[0,1].set_title(r'$\\bf{B}_1$)  Histogram of raw heights')\n",
    "\n",
    "axs[1,1].hist(weight,bins='fd',color=(.5,.5,.5),edgecolor=(.8,.8,.8))\n",
    "axs[1,1].axvline(ind_weight,color='k',linestyle='--',linewidth=2)\n",
    "axs[1,1].set(xlabel='Weight (kg)',ylabel='Count',xlim=[ind_weight-20,ind_weight+20])\n",
    "axs[1,1].set_title(r'$\\bf{B}_2$)  Histogram of raw weights')\n",
    "\n",
    "\n",
    "\n",
    "# plot the z-normalized distributions\n",
    "axs[0,2].hist(height_z,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[0,2].axvline(ind_height_z,color='k',linestyle='--',linewidth=2)\n",
    "axs[0,2].set(xlabel='Normalized height (z)',ylabel='Count',xlim=[-3,3])\n",
    "axs[0,2].set_title(r'$\\bf{C}_1$)  Histogram of z-heights')\n",
    "\n",
    "axs[1,2].hist(weight_z,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[1,2].axvline(ind_weight_z,color='k',linestyle='--',linewidth=2)\n",
    "axs[1,2].set(xlabel='Normalized weight (z)',ylabel='Count',xlim=[-3,3])\n",
    "axs[1,2].set_title(r'$\\bf{C}_2$)  Histogram of z-weights')\n",
    "\n",
    "\n",
    "# export\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_zscore_example.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "LnygNhqQdB30"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "3W-JHyrzlzoK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.2: Z-scoring preserves relative values"
   ],
   "metadata": {
    "id": "klc8JRqw5ki8"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "x = np.random.normal(250,58,size=28)\n",
    "zx = (x-np.mean(x)) / np.std(x,ddof=1)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.plot(x,zx,'ks',markersize=10,markerfacecolor=(.4,.4,.4),alpha=.6)\n",
    "plt.xlabel('Original data')\n",
    "plt.ylabel('Z-transformed data')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_zRelativelyEqual.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "3gcUSz1W5j6C"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "VGJCOsBm5j3a"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.3: Z-scoring a non-normal distribution"
   ],
   "metadata": {
    "id": "fX8XfzOFdBx3"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "N = 5000\n",
    "X = (1+np.random.exponential(size=N))*10\n",
    "Xz = (X-np.mean(X)) / np.std(X,ddof=1)\n",
    "\n",
    "_,axs = plt.subplots(2,1,figsize=(8,7))\n",
    "axs[0].hist(X,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[0].axvline(np.mean(X),color='k',linestyle='--')\n",
    "axs[0].set(xlabel='Raw data value',ylabel='Count',xlim=[8,80])\n",
    "axs[0].set_title(r'$\\bf{A}$)  Histogram of raw data')\n",
    "\n",
    "axs[1].hist(Xz,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[1].axvline(np.mean(Xz),color='k',linestyle='--')\n",
    "axs[1].set(xlabel='Z-transformed data value',ylabel='Count',xlim=[-1.2,5])\n",
    "axs[1].set_title(r'$\\bf{B}$)  Histogram of z-transformed data')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_zscore_positive.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "7HmgMLQcdP3q"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "yn9nLccTTIpm"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.4: Modified-z"
   ],
   "metadata": {
    "id": "QLmGlnPP-vRF"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Some non-normal data\n",
    "x1 = np.random.randn(500) + 2.5\n",
    "x2 = np.random.randn(2500) - 2\n",
    "y  = np.concatenate((x1,x2))\n",
    "y  = y - np.min(y)+3\n",
    "\n",
    "\n",
    "# regular z-score\n",
    "y_z = (y-np.mean(y)) / np.std(y,ddof=1)\n",
    "\n",
    "# modified z\n",
    "MAD = np.median( np.abs(y - np.median(y)) )\n",
    "y_zm = stats.norm.ppf(3/4) * (y - np.median(y)) / MAD\n",
    "\n",
    "\n",
    "# their histograms\n",
    "yy_z,xx_z = np.histogram(y_z,bins='fd')\n",
    "yy_zm,xx_zm = np.histogram(y_zm,bins='fd')\n",
    "\n",
    "\n",
    "\n",
    "_,axs = plt.subplots(1,3,figsize=(10,4))\n",
    "\n",
    "axs[0].hist(y,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[0].set(xlabel='Data value',ylabel='Counts',title=r'$\\bf{A}$)  Original data dist.')\n",
    "\n",
    "axs[1].plot((xx_z[1:]+xx_z[:-1])/2,yy_z,linewidth=2,color='k',label='Regular')\n",
    "axs[1].plot((xx_zm[1:]+xx_zm[:-1])/2,yy_zm,'--',linewidth=2,color=(.4,.4,.4),label='Modified')\n",
    "axs[1].legend(loc='upper right',frameon=False,bbox_to_anchor=[1.08,1.05])\n",
    "axs[1].set(xlabel='Transformed value',ylabel='Counts',title=r'$\\bf{B}$)  Z-score dists.')\n",
    "\n",
    "\n",
    "axs[2].plot(y_z,y_zm,'ko',markerfacecolor='w')\n",
    "xval = np.min( [np.min(y_z),np.min(y_zm)] )\n",
    "yval = np.max( [np.max(y_z),np.max(y_zm)] )\n",
    "axs[2].plot([xval,yval],[xval,yval],'--',color=(.6,.6,.6),label='Unity')\n",
    "axs[2].set(xlabel='\"Regular\" z',ylabel='Modified z',title=r'$\\bf{C}$)  Comparison')\n",
    "axs[2].legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_modVreg_zscore.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "6jgQnLVoRd_U"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "TAFK1iAhfJ3p"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.5: Min-max scaling"
   ],
   "metadata": {
    "id": "BO6e5K5afJ6k"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def minmaxScaling(x):\n",
    "\n",
    "  # compute min/max values\n",
    "  minx = np.min(x)\n",
    "  maxx = np.max(x)\n",
    "\n",
    "  # transformation\n",
    "  xScaled = (x-minx) / (maxx-minx)\n",
    "\n",
    "  # output\n",
    "  return xScaled\n",
    "\n",
    "# Note: I wrote the function over multiple lines for clarity; you could reduce it to one line!"
   ],
   "metadata": {
    "id": "gqUOOSVtr2DO"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "## create some data\n",
    "\n",
    "N = 42\n",
    "data = np.log(np.random.rand(N))*234 + 934\n",
    "\n",
    "# now min-max scale\n",
    "dataS = minmaxScaling(data)\n",
    "\n",
    "\n",
    "# now plot\n",
    "fig,ax = plt.subplots(1,3,figsize=(10,4))\n",
    "randomXoffsets = 1+np.random.randn(N)/20\n",
    "ax[0].plot(randomXoffsets,data,'ks',markerfacecolor='w')\n",
    "ax[0].set(xlim=[0,2],xticks=[],ylabel='Original data scale')\n",
    "ax[0].set_title(r'$\\bf{A}$)  Original data')\n",
    "\n",
    "ax[1].plot(randomXoffsets,dataS,'ks',markerfacecolor='w')\n",
    "ax[1].set(xlim=[0,2],xticks=[],ylabel='Unity-normed data scale')\n",
    "ax[1].set_title(r'$\\bf{B}$)  Scaled data')\n",
    "\n",
    "ax[2].plot(data,dataS,'ks',markerfacecolor='w')\n",
    "ax[2].set(xlabel='Original values',ylabel='Scaled values')\n",
    "ax[2].set_title(r'$\\bf{C}$)  Scatter plot')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_minmax.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "0Y9cQy73fJ9V"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "yP-ulLxdJNch"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.6: Another example of min-max scaling"
   ],
   "metadata": {
    "id": "eZ5tbjLbJNTS"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# generate a Laplace distribution\n",
    "x1 = np.exp(-np.abs(3*np.random.randn(40)))\n",
    "x2 = np.exp(-np.abs(3*np.random.randn(40)))\n",
    "x = (x1-x2)*42 - 13\n",
    "\n",
    "\n",
    "# minmax scale\n",
    "xm = minmaxScaling(x)\n",
    "\n",
    "\n",
    "_,axs = plt.subplots(2,1,figsize=(4,8))\n",
    "\n",
    "axs[0].hist(x,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[0].set(yticks=[],ylabel='Count (a.u.)',xlabel='Raw data values')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Histogram of raw data')\n",
    "\n",
    "axs[1].hist(xm,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[1].set(yticks=[],ylabel='Count (a.u.)',xlabel='Minmax-scaled values')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Hist. of transformed data')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_minmax_exampleHist.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "WOiMYEmXDaAf"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "GfFxBRLo0Upz"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Percent change"
   ],
   "metadata": {
    "id": "AKuw-RgY0UnE"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Note: This code does not correspond to any figure in the book;\n",
    "#  I include it here to illustrate the code.\n",
    "#  Notice in the graph how the 'ref' data value turns into 0.\n",
    "\n",
    "# a range of values for \"new\"\n",
    "new = np.linspace(3,210,31)\n",
    "\n",
    "# reference value\n",
    "ref = 135\n",
    "\n",
    "# compute percent change\n",
    "pctchg = 100*(new-ref) / ref\n",
    "\n",
    "\n",
    "# visualize the transformed data\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.plot(new,pctchg,'ks-',markerfacecolor='m',zorder=30)\n",
    "plt.axhline(0,color='k',linestyle='--')\n",
    "plt.axvline(ref,color='k',linestyle='--')\n",
    "plt.xlabel('Original data values')\n",
    "plt.ylabel('Transformed data values')\n",
    "plt.title('Percent change',loc='center')\n",
    "plt.grid()\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "zUfPl_tv0WER"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "AL4taua6xZvt"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Rank transform"
   ],
   "metadata": {
    "id": "DBrnOFZYxZtJ"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "x = [ 4,1,2,3,3 ]\n",
    "\n",
    "stats.rankdata(x)"
   ],
   "metadata": {
    "id": "_XRPBcd-HpdX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "TpoAsfoJxZqQ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.7: Log transform"
   ],
   "metadata": {
    "id": "hJpZajWYxZlF"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# generate data from power-law distribution\n",
    "X = np.random.normal(0,5,size=2000)**2\n",
    "Xlog = np.log(X)\n",
    "\n",
    "\n",
    "_,axs = plt.subplots(2,2,figsize=(10,7))\n",
    "\n",
    "# scatter plots\n",
    "axs[0,0].plot(X,'ko',markerfacecolor='w')\n",
    "axs[0,0].set(xlabel='Data index',ylabel='Data value')\n",
    "axs[0,0].set_title(r'$\\bf{A}$)  Raw data')\n",
    "\n",
    "axs[0,1].plot(Xlog,'ko',markerfacecolor='w')\n",
    "axs[0,1].set(xlabel='Data index',ylabel='ln(data)')\n",
    "axs[0,1].set_title(r'$\\bf{B}$)  Transformed data')\n",
    "\n",
    "\n",
    "# histograms\n",
    "axs[1,0].hist(X,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[1,0].set_title(r'$\\bf{C}$)  Histogram of raw data')\n",
    "\n",
    "axs[1,1].hist(Xlog,bins='fd',color=(.9,.9,.9),edgecolor='k')\n",
    "axs[1,1].set_title(r'$\\bf{D}$)  Histogram of ln(data)')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_logdata_example.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "56LRT3ClRwgw"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "gB9xF470fKAF"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.8: Logarithm vs square root"
   ],
   "metadata": {
    "id": "FZSXe1UVfKDE"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# some right-tailed non-normally distributed data\n",
    "X = np.random.randn(10000)**2\n",
    "\n",
    "y_r,x_r = np.histogram(X,bins=100) # r = raw\n",
    "y_l,x_l = np.histogram(np.log(X),bins=100) # l = log\n",
    "y_s,x_s = np.histogram(np.sqrt(X),bins=100) # s = sqrt\n",
    "\n",
    "\n",
    "\n",
    "_,axs = plt.subplots(1,3,figsize=(10,4))\n",
    "\n",
    "# theory\n",
    "q = np.linspace(.1,10,100)\n",
    "axs[0].plot(q,np.log(q),'k',linewidth=3,label='ln(x)')\n",
    "axs[0].plot(q,np.sqrt(q),'--',color='gray',linewidth=3,label='sqrt(x)')\n",
    "axs[0].set(xlim=[0,10],xlabel='Raw data value',ylabel='Transformed data value')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Transformation')\n",
    "axs[0].legend()\n",
    "\n",
    "# untransformed data\n",
    "axs[1].plot((x_r[1:]+x_r[:-1])/2,y_r,'k',linewidth=2)\n",
    "axs[1].set(xlabel='Raw data value',ylabel='Count')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Data histogram')\n",
    "\n",
    "axs[2].plot((x_l[1:]+x_l[:-1])/2,y_l,'k',linewidth=2,label='ln(x)')\n",
    "axs[2].plot((x_s[1:]+x_s[:-1])/2,y_s,'--',color='gray',linewidth=2,label='sqrt(x)')\n",
    "axs[2].set(xlabel='Transformed data value',ylabel='Count')\n",
    "axs[2].set_title(r'$\\bf{C}$)  Trans data hist.')\n",
    "axs[2].legend()\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_logsqrt.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "VrWgPt92Tngk"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "naDAB9CV_YWa"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.9: Fisher-z transform"
   ],
   "metadata": {
    "id": "1A3o_eQ7_YTs"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "N = 1000\n",
    "\n",
    "# uniform data in range [-1,1]\n",
    "Y = np.random.uniform(-1,1,size=N)\n",
    "# Note: Fisher-z is undefined for Y==|1|, but the probability of that happening in\n",
    "#       random uniform data is so vanishingly small that coding an exception is unnecessary.\n",
    "\n",
    "# transform\n",
    "fY = np.arctanh(Y)\n",
    "\n",
    "_,axs = plt.subplots(1,3,figsize=(10,np.pi))\n",
    "\n",
    "axs[0].hist(Y,30,color=(.8,.8,.8),edgecolor='k')\n",
    "axs[0].set(xlabel='Data values', ylabel='Count')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Raw data hist.')\n",
    "\n",
    "axs[1].hist(fY,30,color=(.8,.8,.8),edgecolor='k')\n",
    "axs[1].set(xlim=[-5,5],xlabel='Data values', ylabel='Count')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Fisher data hist.')\n",
    "\n",
    "axs[2].plot(Y,fY,'k.')\n",
    "axs[2].set_title(r'$\\bf{C}$)  Transformation')\n",
    "axs[2].set(xlabel='Original data',ylabel='Transformed data')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_fisherz.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "BFYuHnBb_aJN"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "XJ0DKAtsRTtK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.10: Transform any distribution to Gaussian"
   ],
   "metadata": {
    "id": "D0Ib-xIgRTqA"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "X = np.random.normal(0,5,size=2000)**2\n",
    "\n",
    "X_r = stats.rankdata( X )\n",
    "a,b = -.999,.999\n",
    "X_r2 = (X_r-np.min(X_r)) / (np.max(X_r)-np.min(X_r)) * (b-a) + a\n",
    "X_t = np.arctanh( X_r2 )\n",
    "\n",
    "_,axs = plt.subplots(1,3,figsize=(10,4))\n",
    "axs[0].hist(X,40,color=(.8,.8,.8),edgecolor='k')\n",
    "axs[0].set(xlabel='Data value',ylabel='Count')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Raw data hist.')\n",
    "\n",
    "axs[1].hist(X_t,40,color=(.8,.8,.8),edgecolor='k')\n",
    "axs[1].set(xlabel='Data value',ylabel='Count')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Trans. data hist.')\n",
    "\n",
    "axs[2].plot(X,X_t,'ko',markersize=8,markerfacecolor=(.8,.8,.8))\n",
    "axs[2].set(xlabel='Original data',ylabel='Transformed data')\n",
    "axs[2].set_title(r'$\\bf{C}$)  Data comparison')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_any2gauss.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "BOvTlG_ZRWQJ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "FAzaJE4z_YOJ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 6.11: Linearization of nonlinear relationship"
   ],
   "metadata": {
    "id": "Lk3dSt460-6r"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# sample size\n",
    "N = 200\n",
    "\n",
    "# the data (note how the nonlinearity is implemented)\n",
    "x = np.linspace(.01,3,N)\n",
    "y = np.exp(x) + np.random.randn(N)*np.linspace(.1,1,N)*3\n",
    "y = np.abs(y)\n",
    "\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "axs[0].plot(x,y,'ko',markerfacecolor='w')\n",
    "axs[0].set(xlabel='x',ylabel='y')\n",
    "axs[0].set_title(r'$\\bf{A})$  Original data')\n",
    "\n",
    "axs[1].plot(x,np.log(y),'ko',markerfacecolor='w')\n",
    "axs[1].set(xlabel='x',ylabel='ln(y)')\n",
    "axs[1].set_title(r'$\\bf{B})$  Transformed data')\n",
    "\n",
    "a,b = np.polyfit(x,np.log(y),1)\n",
    "axs[1].plot(x,a*x+b,'k',linewidth=3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_linearizedFit.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "QgRVs54w1T8S"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "t-C3sLlDTnd9"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 1"
   ],
   "metadata": {
    "id": "fCCvq5xTTnbH"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# a function to scale to any arbitrary bounds\n",
    "def my_minmaxScalar(x,newmin,newmax):\n",
    "\n",
    "  # get min and max\n",
    "  minx = np.min(x)\n",
    "  maxx = np.max(x)\n",
    "\n",
    "  # (intermediate) min-max scale to [0,1]\n",
    "  xS = (x-minx) / (maxx-minx)\n",
    "\n",
    "  # scale to [a,b]\n",
    "  xSS = xS*(newmax-newmin) + newmin\n",
    "\n",
    "  return xSS\n",
    "\n",
    "  # this solution harder to read, but is what your equation should look like on paper:\n",
    "  # return (x-np.min(x)) / (np.max(x)-np.min(x)) * (newmax-newmin) + newmin"
   ],
   "metadata": {
    "id": "eXro4KyI-HoK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# test it on some data\n",
    "data = np.random.randn(10)\n",
    "a,b = 14.3,34\n",
    "\n",
    "# apply the transform\n",
    "y = my_minmaxScalar(data,a,b)\n",
    "\n",
    "\n",
    "print(f'Min value: {np.min(y)}')\n",
    "print(f'Max value: {np.max(y)}')"
   ],
   "metadata": {
    "id": "TtbPWcjs7cZP"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# import scalar class from sklearn\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "# define an instance of this class, with boundary parameters\n",
    "scaler = MinMaxScaler(feature_range=(a,b))\n",
    "\n",
    "# fit to the data (but the data need to be a column vector)\n",
    "scaler.fit(data.reshape(-1,1))\n",
    "\n",
    "# finally, call the transformer object with the data (again as column vector)\n",
    "y2 = scaler.transform(data.reshape(-1,1))\n",
    "\n",
    "# print the min and max values\n",
    "print(f'Min value: {np.min(y2)}')\n",
    "print(f'Max value: {np.max(y2)}')"
   ],
   "metadata": {
    "id": "HQo57j0dwOjd"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# print out all individual numbers to confirm they're the same\n",
    "np.hstack((y.reshape(-1,1),y2.reshape(-1,1)))"
   ],
   "metadata": {
    "id": "PL9XZTnMK20q"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# scatter plot\n",
    "plt.plot(y,y2,'ko',markerfacecolor=(.1,.2,.3),markersize=12)\n",
    "plt.xlabel('My scalar')\n",
    "plt.ylabel(\"sklearn's scalar\")\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "IK6RVQVewOl_"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "ZbJOcN-gOS6-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 2"
   ],
   "metadata": {
    "id": "8AHdKDCXOS32"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# this is the only line that's modified from the code earlier in the chapter\n",
    "X = (np.random.randn(10000)+3)**2\n",
    "\n",
    "\n",
    "\n",
    "y_r,x_r = np.histogram(X,bins=100) # r = raw\n",
    "y_l,x_l = np.histogram(np.log(X),bins=100) # l = log\n",
    "y_s,x_s = np.histogram(np.sqrt(X),bins=100) # s = sqrt\n",
    "\n",
    "\n",
    "\n",
    "_,axs = plt.subplots(1,3,figsize=(10,4))\n",
    "\n",
    "# theory\n",
    "q = np.linspace(.1,10,100)\n",
    "axs[0].plot(q,np.log(q),'k',linewidth=3,label='ln(x)')\n",
    "axs[0].plot(q,np.sqrt(q),'--',color='gray',linewidth=3,label='sqrt(x)')\n",
    "axs[0].set(xlim=[0,10],xlabel='Raw data value',ylabel='Transformed data value')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Transformation')\n",
    "axs[0].legend()\n",
    "\n",
    "# untransformed data\n",
    "axs[1].plot((x_r[1:]+x_r[:-1])/2,y_r,'k',linewidth=2)\n",
    "axs[1].set(xlabel='Raw data value',ylabel='Count')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Data histogram')\n",
    "\n",
    "axs[2].plot((x_l[1:]+x_l[:-1])/2,y_l,'k',linewidth=2,label='ln(x)')\n",
    "axs[2].plot((x_s[1:]+x_s[:-1])/2,y_s,'--',color='gray',linewidth=2,label='sqrt(x)')\n",
    "axs[2].set(xlabel='Transformed data value',ylabel='Count')\n",
    "axs[2].set_title(r'$\\bf{C}$)  Trans data hist.')\n",
    "axs[2].legend()\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "Z6FyGBLpOUWI"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# data\n",
    "X1 = (np.random.randn(1000)+0)**2\n",
    "X2 = (np.random.randn(1000)+3)**2\n",
    "\n",
    "# sqrt transform\n",
    "X1s = np.sqrt(X1)\n",
    "X2s = np.sqrt(X2)\n",
    "\n",
    "\n",
    "_,axs = plt.subplots(2,2,figsize=(10,8))\n",
    "\n",
    "# QQ plots\n",
    "stats.probplot(X1,plot=axs[0,0],fit=True)\n",
    "stats.probplot(X2,plot=axs[0,1],fit=True)\n",
    "\n",
    "stats.probplot(X1s,plot=axs[1,0],fit=True)\n",
    "stats.probplot(X2s,plot=axs[1,1],fit=True)\n",
    "\n",
    "for a in axs.flatten():\n",
    "  a.get_lines()[0].set_markerfacecolor('k')\n",
    "  a.get_lines()[0].set_markeredgecolor('k')\n",
    "  a.get_lines()[1].set_color('gray')\n",
    "  a.get_lines()[1].set_linewidth(3)\n",
    "  a.set_title(' ')\n",
    "  a.set_ylabel('Data values (sorted)')\n",
    "\n",
    "\n",
    "axs[0,0].set_title(r'$\\bf{A}$)  QQ plot of $x^2$')\n",
    "axs[0,1].set_title(r'$\\bf{B}$)  QQ plot of $\\sqrt{x^2}$')\n",
    "axs[1,0].set_title(r'$\\bf{C}$)  QQ plot of $(x+3)^2$')\n",
    "axs[1,1].set_title(r'$\\bf{D}$)  QQ plot of $\\sqrt{(x+3)^2}$')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_ex2.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "idmG-YOICkUc"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# repeat for log transform\n",
    "X1s = np.log(X1)\n",
    "X2s = np.log(X2)\n",
    "\n",
    "\n",
    "_,axs = plt.subplots(2,2,figsize=(10,8))\n",
    "\n",
    "# QQ plots\n",
    "stats.probplot(X1,plot=axs[0,0],fit=True)\n",
    "stats.probplot(X2,plot=axs[0,1],fit=True)\n",
    "\n",
    "stats.probplot(X1s,plot=axs[1,0],fit=True)\n",
    "stats.probplot(X2s,plot=axs[1,1],fit=True)\n",
    "\n",
    "for a in axs.flatten():\n",
    "  a.get_lines()[0].set_markerfacecolor('k')\n",
    "  a.get_lines()[0].set_markeredgecolor('k')\n",
    "  a.get_lines()[1].set_color('gray')\n",
    "  a.get_lines()[1].set_linewidth(3)\n",
    "  a.set_title(' ')\n",
    "  a.set_ylabel('Data values (sorted)')\n",
    "\n",
    "\n",
    "axs[0,0].set_title(r'$\\bf{A}$)  QQ plot of $x^2$')\n",
    "axs[0,1].set_title(r'$\\bf{B}$)  QQ plot of $\\ln{x^2}$')\n",
    "axs[1,0].set_title(r'$\\bf{C}$)  QQ plot of $(x+3)^2$')\n",
    "axs[1,1].set_title(r'$\\bf{D}$)  QQ plot of $\\ln{(x+3)^2}$')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "y-uhGzhdCkR1"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Shift the data to be mean-centered\n",
    "X = (np.random.randn(10000)+3)**2\n",
    "Y = np.sqrt(X)\n",
    "Y -= np.mean(Y)\n",
    "y_s2,x_s2 = np.histogram(Y,bins=100)\n",
    "\n",
    "plt.plot((x_s[1:]+x_s[:-1])/2,y_s,linewidth=2,label='sqrt(x)')\n",
    "plt.plot((x_s2[1:]+x_s2[:-1])/2,y_s2,linewidth=2,label='shifted sqrt(x)')\n",
    "plt.xlabel('Transformed data value')\n",
    "plt.ylabel('Count')\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "ATq3eCb5OS0Q"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "2ZE_VKPvwOoP"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 3"
   ],
   "metadata": {
    "id": "R4MD7E6n4VXh"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# After a bit of algebra, you should arive at the equation:\n",
    "#  zs+m = x\n",
    "# where: z is the z-transformed data\n",
    "#        s is the standard deviation of the original data\n",
    "#        m is the mean (mu) of the original data\n",
    "#        x is the reconstructed data\n",
    "#\n",
    "# The implication is Yes, the z-score transformation is invertible,\n",
    "#    but only if you store the mean and std of the pre-transformed data."
   ],
   "metadata": {
    "id": "b3ccPdPE4VaO"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "### Now for the empirical illustration:\n",
    "\n",
    "# generate data\n",
    "data = np.random.randint(4,15,25)\n",
    "\n",
    "# pre-transformed descriptive statistics\n",
    "orig_m = np.mean(data)\n",
    "orig_s = np.std(data,ddof=1)\n",
    "\n",
    "# transform to z-score\n",
    "dataz = (data-orig_m) / orig_s\n",
    "\n",
    "# back-transform\n",
    "data_inv = dataz*orig_s + orig_m\n",
    "\n",
    "\n",
    "# There are several ways to confirm the equality of data and data_inv.\n",
    "#   - You can print both vectors and inspect them.\n",
    "#   - You can subtract them to see that their difference is zero.\n",
    "#   - You can plot them on top of each other and see that they overlap (what I do below).\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(data,'ks',markersize=12,markerfacecolor='w',label='Original')\n",
    "plt.plot(data,'kx',markersize=5,markerfacecolor='w',label='Back-transformed')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel('Data index')\n",
    "plt.ylabel('Data value')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "hzGkN8Dt3tbW"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "ppsYQf3SxmtO"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 4"
   ],
   "metadata": {
    "id": "sVC04GFl_dSr"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# parameters\n",
    "N = 313\n",
    "lo_bnd = 3*np.pi\n",
    "hi_bnd = np.exp(np.pi)\n",
    "\n",
    "# create data\n",
    "data = np.random.uniform(lo_bnd,hi_bnd,size=N)\n",
    "\n",
    "### my algorithm\n",
    "# step 1: scale to [-1,1] (but not exactly -1 or +1 b/c of domain restrictions)\n",
    "dataScale = my_minmaxScalar(data,-.999,.999) # using function from ex1\n",
    "\n",
    "# step 2: Fisher-z to transform to normal\n",
    "dataFish = np.arctanh(dataScale)\n",
    "\n",
    "# step 3: shift mean\n",
    "# (min/max scaling doesn't guarantee a particular mean,\n",
    "#  so I mean-centered the data before shifting to the original mean)\n",
    "dataFinal = dataFish - np.mean(dataFish) + (lo_bnd+hi_bnd)/2\n",
    "\n",
    "# report the expected and empirical means\n",
    "print(f'Expected mean : {(lo_bnd+hi_bnd)/2:.3f}')\n",
    "print(f'Empirical mean: {np.mean(dataFinal):.3f}')"
   ],
   "metadata": {
    "id": "M2wvHmv7CDv4"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# create histogram and plot the distributions\n",
    "\n",
    "# distribution values\n",
    "yO,xO = np.histogram(data,bins='fd')\n",
    "yF,xF = np.histogram(dataFinal,bins='fd')\n",
    "\n",
    "# bin centers\n",
    "xO = (xO[1:]+xO[:-1])/2\n",
    "xF = (xF[1:]+xF[:-1])/2\n",
    "\n",
    "# plot the histograms\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(xO,yO,'ks--',linewidth=3,label='Original data')\n",
    "plt.plot(xF,yF,'o-',color=(.7,.7,.7),linewidth=3,label='Transformed data')\n",
    "\n",
    "# overlay the means\n",
    "plt.axvline(x=np.mean(data),color=(.8,.8,.8), linestyle=':',label='Average')\n",
    "\n",
    "# flotsam and jetsam\n",
    "plt.xlabel('Data value')\n",
    "plt.ylabel('Count')\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_ex4.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "vDk4EdiyANEh"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Observation: It's interesting that the Fisher-z transform often stretches out\n",
    "# the data compared to its original uniform distribution, but here the transformed\n",
    "# distribution is narrower. That happens because the data are first compressed to a\n",
    "# range of [-1,1], and most of the data values lie towards zero.\n",
    "#\n",
    "# Oftentimes in statistics, the shape of the distribution is more important than\n",
    "# the numerical range of the distribution."
   ],
   "metadata": {
    "id": "_2J-0IkfCD1e"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "n_NbTcmuCD4a"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 5"
   ],
   "metadata": {
    "id": "V6S6hVC9CD7S"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# z-score a simple array\n",
    "data = np.arange(3,10)\n",
    "\n",
    "# compute z-score manually and via scipy.stats\n",
    "data_z_me = (data-np.mean(data)) / np.std(data,ddof=1)\n",
    "data_z_sp = stats.zscore(data,ddof=1) # note the ddof parameter!\n",
    "\n",
    "# print the results to confirm they match\n",
    "print(data_z_me)\n",
    "print(data_z_sp)"
   ],
   "metadata": {
    "id": "mgpv2r0Tsgsj"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# z-score columns vs whole matrix.\n",
    "\n",
    "data = np.arange(0,12).reshape(4,3)**2\n",
    "\n",
    "data_z_col = stats.zscore(data,ddof=1,axis=0) # axis=0 is the default; I'm writing here so it's clear to you\n",
    "data_z_mat = stats.zscore(data,ddof=1,axis=None) # gets the mean/std from the whole matrix\n",
    "\n",
    "print('Original data matrix:')\n",
    "print(data)\n",
    "\n",
    "print(' ')\n",
    "print('Column-wise z-scoring:')\n",
    "print(data_z_col)\n",
    "\n",
    "print(' ')\n",
    "print('Matrix-wise z-scoring:')\n",
    "print(data_z_mat)"
   ],
   "metadata": {
    "id": "wrHVhRmnsguu"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "GVgPKnrCCEBB"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 6"
   ],
   "metadata": {
    "id": "IwEvt1xqCED4"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# The insight here is that the Fisher-z transform is nonlinear,\n",
    "# and has a bigger impact on data that are closer to |1|.\n",
    "# This means you can create skew from a uniform distribution\n",
    "# with asymmetric boundaries."
   ],
   "metadata": {
    "id": "PDEquB11wxEr"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# no skew\n",
    "X1 = np.arctanh(np.random.uniform(-.999,.999,size=5000))\n",
    "\n",
    "# strong positive skew\n",
    "X2 = np.arctanh(np.random.uniform(-.2,.999,size=5000))\n",
    "\n",
    "# slight negative skew\n",
    "X3 = np.arctanh(np.random.uniform(-.999,.8,size=5000))\n",
    "\n",
    "\n",
    "# compute the empirical skews\n",
    "skews = [0]*3\n",
    "skews[0] = stats.skew(X1)\n",
    "skews[1] = stats.skew(X2)\n",
    "skews[2] = stats.skew(X3)\n",
    "\n",
    "\n",
    "_,axs = plt.subplots(1,3,figsize=(10,3.5))\n",
    "\n",
    "axs[0].hist(X1,40,color=(.8,.8,.8),edgecolor='k')\n",
    "axs[0].set_title(f'$\\\\bf{{A}}$)  skew = {skews[0]:.2f}')\n",
    "\n",
    "axs[1].hist(X2,40,color=(.8,.8,.8),edgecolor='k')\n",
    "axs[1].set_title(f'$\\\\bf{{B}}$)  skew = {skews[1]:.2f}')\n",
    "\n",
    "axs[2].hist(X3,40,color=(.8,.8,.8),edgecolor='k')\n",
    "axs[2].set_title(f'$\\\\bf{{C}}$)  skew = {skews[2]:.2f}')\n",
    "\n",
    "\n",
    "for a in axs:\n",
    "  a.set(yticks=[],ylabel='Counts',xticks=[],xlabel='Data values')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_ex6.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "kUmUQUuPsoYX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "sTX6_oE0soeQ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 7"
   ],
   "metadata": {
    "id": "2eza1TxJsog_"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# generate the data\n",
    "sigmas = np.linspace(.1,1,20)\n",
    "X = np.random.randn(13524)\n",
    "\n",
    "# initialize the data results matrices\n",
    "M = np.zeros((len(sigmas),2))\n",
    "data = [0]*len(sigmas)\n",
    "\n",
    "# compute and store all moments in a matrix\n",
    "for i,s in enumerate(sigmas):\n",
    "\n",
    "  # create the data\n",
    "  data[i] = np.exp(X*s)\n",
    "\n",
    "  # compute its stadard deviation\n",
    "  M[i,0] = np.std( data[i] ,ddof=1)\n",
    "\n",
    "  # compute its MAD\n",
    "  M[i,1] = np.median( np.abs(data[i] - np.median(data[i])) )"
   ],
   "metadata": {
    "id": "PmicKMdtsoj3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "_,axs = plt.subplots(1,2,figsize=(10,5))\n",
    "\n",
    "# plot the moments\n",
    "m = 'osp*' # markers\n",
    "l = ['-','--',':','-.'] # line type\n",
    "c = [0,.2,.4,.6] # grayscale\n",
    "\n",
    "\n",
    "# plot the std and MAD (plotted first for the vertical line heights)\n",
    "axs[1].plot(sigmas,M[:,0],'ks-',markersize=10,markerfacecolor='w',label='std')\n",
    "axs[1].plot(sigmas,M[:,1],'ko-',markersize=10,markerfacecolor=(.8,.8,.8),label='MAD')\n",
    "axs[1].plot(sigmas,M[:,1]/stats.norm.ppf(3/4),'kx-',markersize=10,markerfacecolor=(.8,.8,.8),label='Scaled MAD')\n",
    "axs[1].legend()\n",
    "axs[1].set(xlabel=r'$\\sigma$',ylabel='Dispersion value')\n",
    "axs[1].set_title(r'$\\bf{B)}$  Standard deviation and MAD')\n",
    "\n",
    "\n",
    "\n",
    "# now plot selected distributions\n",
    "for i in np.linspace(0,len(sigmas)-1,5).astype(int):\n",
    "\n",
    "  # get histogram values\n",
    "  y,x = np.histogram(data[i],bins='fd')\n",
    "\n",
    "  # plot as line\n",
    "  h = axs[0].plot((x[:-1]+x[1:])/2,y,'.-',linewidth=2,markersize=8,\n",
    "              label=f'$\\sigma={sigmas[i]:.2f}$')\n",
    "\n",
    "  # add the vertical lines in the moments plot\n",
    "  axs[1].annotate(f'{sigmas[i]:.2f}',xy=(sigmas[i],.1),horizontalalignment='center',fontsize=10,\n",
    "                  xytext=(sigmas[i],axs[1].get_ylim()[1]/2),arrowprops=dict(width=4,linewidth=0,color=h[0].get_color(),alpha=.8))\n",
    "\n",
    "\n",
    "# some niceties\n",
    "axs[0].legend()\n",
    "axs[0].set(xlim=[0,6],xlabel='Data value',ylabel='Count')\n",
    "axs[0].set_title(r'$\\bf{A)}$  Distributions')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_ex7.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "uT1AeS0ksomj"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "JfKKp1FXso5-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 8"
   ],
   "metadata": {
    "id": "SayqLaJ1KMHt"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create two non-normal distributions\n",
    "N = 300\n",
    "sample1 = np.random.randn(N)**2\n",
    "sample2 = np.random.randn(N)**2\n",
    "\n",
    "# and their difference\n",
    "difference = sample1 - sample2\n",
    "\n",
    "# compute their histograms\n",
    "y1,x1 = np.histogram(sample1,bins='fd')\n",
    "y2,x2 = np.histogram(sample2,bins='fd')\n",
    "yd,xd = np.histogram(difference,bins='fd')\n",
    "\n",
    "# and plot\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "axs[0].plot((x1[:-1]+x1[1:])/2,y1,'k',linewidth=2,label='Sample 1')\n",
    "axs[0].plot((x2[:-1]+x2[1:])/2,y2,'--',color='gray',linewidth=2,label='Sample 2')\n",
    "axs[0].set(xlabel='Data value',ylabel='Count')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Individual histograms')\n",
    "axs[0].legend()\n",
    "\n",
    "axs[1].plot((xd[:-1]+xd[1:])/2,yd,'k',linewidth=2)\n",
    "axs[1].set(xlabel='Difference value',ylabel='Count')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Difference histograms')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('trans_ex8.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "8pBE3UawKNTS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "gcU5UydJKMCr"
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}